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Abstract: We perform dynamical QCD simulations with nj = 2 overlap fermions by 
hybrid Monte-Carlo method on 6^ to 8^ x 16 lattices. We study the problem of topolog- 
ical sector changing. A new method is proposed which works without topological sector 
changes. We use this new method to determine the topological susceptibility at various 
quark masses. 
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1. Introduction 

The overlap operator [1,2] , gives a theoretically sound solution of the chirality problem on 
the lattice. It satisfies the Ginsparg- Wilson relation [3,4], which ensures the exact chiral 
symmetry at finite lattice spacing [5], moreover the difi'erence in the number of left and 
right handed zero modes can be taken as a definition of the topological charge {Q) which 
gives the correct result in the continuum limit [6]. 

However, the numerical implementation of dynamical overlap fermions is still a great 
challenge today (for early studies with dynamical overlap fermions we refer to [7-9]). The 
presence of nested conjugate gradients for the inversion of the Dirac operator makes the 
simulations considerably slower than simulations with Wilson fermions. Furthermore one 
has to face the non-continuity of the fermion determinant at the boundary of topological 
sectors. This additional difficulty can be treated exactly in frame of the Hybrid Monte Carlo 
(HMC) algorithm by modifying the molecular dynamics trajectory at the boundary [10]. 
Clearly the crossing rate between different topological sectors is heavily affected by this 
modification. Inappropriate treatment might confine the system into a certain topological 
sector which yields an unacceptably large autocorrelation time for Q in the simulation. 

A few exploratory studies are already available in QCD with dynamical overlap fer- 
mions [10-15]. All handle the modification of the trajectory at the boundary in a similar 
style. The original proposal of [10] is modified in [13] in such a way that the acceptance 
rate is increased. It is shown in [14], that the introduction of several pseudofermion fields 
which approximate the fermion determinant, can enhance the crossing rate. The relation 
between the pseudofermionic (over)estimation and rare topological sector changes^ was 
pointed out in [15]. 

^In the staggered formulation there was already a concern that the pseudofermion estimator obstruct 
the change of topological charge [16]. 
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In this paper wc study the problem of changing topological sectors in the case of the 
overlap operator with n/ = 2 fermions in dynamical HMC simulations. Sec. 2 will give a 
short introduction to the sector changing problem for overlap fermions (by summarizing and 
extending the work of [15]), and answers the question why the present treatment is unlikely 
to change Q. In Sec. 3 a new measurement method of expectation values is proposed, 
which circumvents the crossing problem entirely by making simulations constrained to 
fixed topological sectors. In Sec. 4 we present numerical results using the new measurement 
method. In Sec. 5 conclusions are given. 

2. Topological sector changing problem 

After introducing pseudofermion fields [17], our partition function reads: 

Z = j [d?7]e-^« det H'^ = j W] [#^] [#]e-^«-'^^-^"''^, (2.1) 
where H is the hermitian overlap operator: 

with Hq being the massless hermitian overlap operator: 

i^o = mo[75 + sgn(ifw^)]. (2.3) 

Here H\/\r is the standard Wilson operator with negative mass —mo. The sgn function in 
Hq causes a T)irac-6 type singularity in the equation of motion of the momenta of the link 
variables. The (5-f unction gives a contribution whenever an eigenvalue A of Hw changes 
sign. This subspace of the configuration space coincide with topological sector boundaries. 
The reason for this is that in the case of the overlap operator the topological charge is: 

Q= 2^TV(/fo) = ^l]sgn(A,) (2.4) 

i 

This means, that there are potential walls, non-differentiable steps in the action at the 
topological sector boundaries. The reflection-refraction method suggested in [10] handles 

these potential walls correctly. Let's denote the momenta by p and the normal vector of 
the topological sector boundary by n. According to this method one has to modify the 
momenta, when arriving at a potential wall: 

P^Ip- ^('^^p) + if i^^P)^ > (refraction) 

[p- 2n{n,p), if (n,p)^ < 2AS (reflection) 

Thus the trajectory will go through the topological sector boundary only if {n,p)^ > 2 AS. 
In a HMC algorithm {n,p) = 0(1) and has exponential distribution. AS*, however, is not 
the exact value of the height of the potential wall, but it is the change of the pseudofermionic 
action at the boundary. From now on, we will distinguish between these two quantities. 
We call the former A6'exact and the latter A6'pf. 
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Let us take a closer look^ on the relation between AS'cxact and (A5pf). In particular, 
we show that the jump in the pseudofcrmionic action overestimates A5cxact- Let us assume 
that the trajectory crosses the boundary. Let i7_ and iJ+ be the overlap operator evaluated 
on the two sides of the boundary right before and after the crossing, respectively. Clearly 
i7_ and contain the same gauge configuration, but they differ, since one eigenvalue of 
Hw changes sign on the boundary. In the HMC algorithm one chooses the pseudofermion 
field as 

where 77, rj^ are random vectors with Gaussian distribution, in order to generate (/), (p^ with 
the correct distribution. (In a real simulation one chooses new pseudofermion configurations 
only at the beginning of each trajectory, but for simplicity let's consider, that (f) and (f)^ are 
refreshed when hitting the boundary.) The jump of the pseudofermionic action now reads: 

The relation between A^exact and A5pf can be obtained by the following straightforward 
calculation: 



e 



-AiSexact 



detH^ /[dr?t][(i,y]e-'7t'? 

The inequality in the second line is a consequence of the concavity of the function. So 
we conclude to: 

(AS'pf) > AS'exact- 

We can examine this relation in realistic simulations, if we take into account, that 
there is a simple relation between and H^. Let's denote by Aq the eigenvalue of 
which crosses zero at the boundary, and by |0) the eigenvector belonging to Aq. With this 
notation: 

/f+ = iJ_ + c|0)(0|, 

where 

TTl 

c = AsgnAo mo{l 



2mo 

with AsgnAo = ±2 being the jump of sgnAo on the boundary. The expectation value of 
the discontinuity in the pseudofermionic action is: 

(A5pf) = {r^HH.HfH^ - l)r/)^t, = Tt{H^H^^H^ - 1) = 

= TV((1 - c\OmH-'){l - c H~'\0){0\) - 1) = -2cm-'\0) + c'{0\Hf\0). (2.6) 

In a similar way one can get a simple formula for the exact value of the jump on the 
boundary: 

g-ASexact ^ dctff^ ^ 1 ^ 1 ^ 1 ^2 7) 

dct//2 dct(i7+ii/_)2 det(l-ci7+i|0)(0|)2 (1 ~ c(0|F7^|0))2 ' ^ ' 



^The following considerations in this section have already partially appeared in [15]. 
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eq. ( |2.6D and eq. (|2.7D offers a numerically fast way to determine both action jumps, since 
one needs only one inversion of the overlap operator to obtain both of them. 




Figure 1: The jump in the exact vs. pseudofermionic action at /3 = 4.05 and m = 0.1, 0.2. Since the 
average of {n,p)'^ is around ~ 1, topological sector changing would happen considerably frequently 
using ^oxact, than with Sp{. We also indicated the probability of topological sector changing with 
the pseudofermionic action, and an estimate on the probability using the exact action (assuming 
that the two algorithms would behave the same way except for the boundaries). 

For illustration we made a scatter plot (Fig. |l|) from a 6'^ lattice at two different 
masses. (Details of our action will be described in Section ^) One can clearly see, that 
the use of the pseudofermions has an awkward consequence: there are a huge amount 
of crossings, where the topological sector changing fails only due to the overestimation. 
One way to cure this is to use several pseudofermion estimators instead of one [14] . More 
pseudofermions mean smaller spread of the pseudofermionic action distribution, therefore 
the overestimation is smaller, too. However the computational time also increases with 
the number of extra fields. Obviously the best would be to use the exact action in the 
simulations, but only its discontinuity on the boundary can be calculated easily. The next 
section will present a technique, which uses this discontinuity to get the relative weight of 
topological sectors. 

3. The new method 

In this section we propose a new method for the calculation of physical observables by which 
it becomes possible to circumvent the problem of topological sector changing described in 
the previous section. Let us write the partition function in the form (assuming a vanishing 
6 parameter): 
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oo 

Q = — OD 

where Zq is the partition function of the topological sector Q. The expectation value of 
an observable: 

where the restricted expectation value {0)q is 

{0)q = [dU]QO[U]detH^exp{-Sg). 

For reasons which will be clear later the integration goes not only over the configurations 
with Q charge, but also over the boundary of the topological sector as well (though the 
boundary has only zero measure in this case). When calculating the partition function 
in a given topological sector the following boundary prescription is used: we define the 
determinant on the boundary as the limit of determinants approaching the wall from the 
Q side (det Hq). If the measurement of the quantities Zq^i/Zq would be possible, then 
we could recover Zq/Zq for any Q. With these in hand, we would need only the restricted 
expectation values {0)q, whose measurement doesn't require topological sector changings. 

Now we will show a way to measure Zq^i/Zq. It will make use of the fact, that 
we can calculate easily ASgxact on the boundary of topological sectors (see eq. ( |2.7] )). 
The pseudofermionic action is only used to generate configurations in fixed topological 
sectors, so its bad distribution for the jump of the action will not effect us. (In the 
following formulae AS" will automatically mean A5exact-) The main idea is the following: 
an observable measured in sector Q is inversely proportional to Zq and an observable in 
Q + 1 is to Zq^i. If the observables in the two sectors are concentrated only to the common 
wall separating the two sectors, then from the ratio of the two expectation values one can 
recover the ratio of the two sectors. 

First let us measure in the Q sector an operator, which is concentrated to the boundary: 

{dQ,Q+iF)Q = [dU]Q6Q,Q+iF[U]detHl^expi-Sg), (3.1) 

where we introduced the distribution 6q^q^i, a Dirac-(5, which is equal to zero everywhere 
but on the Q,Q + 1 boundary. Then let us measure another operator G on the same wall 
(thus on the boundary separating sectors Q and Q + 1), but now from the Q + 1 sector: 

(5q,q+iG)q+i = / [dU]Q+i5Q,Q+iG[U] det i?^ exp(-5,). (3.2) 

Zq+1 J 

The wall is the same (i.e. [dU]Q5Q^Q+i = [df/]Q+i(^Q,(3+i) in both cases, however due to our 
boundary prescription the determinants are different on it. Therefore if F and G satisfies: 

F[U]detHl,[U] = G[U] det H^Q^,[U], (3.3) 
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then the ratio of eq. (3A) and eq. (|3.2[) gives us 

{^Q,Q+iF)q ^ Zq^i 
{5q,q+iG)q+\ Zq 



(3.4) 



The easiest choice is G{U) = 1 and F{U) = det Hq_^_^/ det Hq = exp(— A5), the ratio 
of sectors becomes: 

y /y (^Q,Q+l eXp(-A5'))Q 

Zq+i/Zq = ^ . (3.5) 

This choice is still not optimal, since the measurement of the numerator is problematic, if 
the distribution of A5 extends to negative values. The exponential function amplifies the 
small fluctuations in the negative AS region, which can destroy the whole measurement: 
a very small fraction of the configurations will dominate the result. As a consequence 
one ends up with relatively large statistical uncertainties. With a slightly different choice 
of F and G we can improve on the situation. With F(U) = @{AS — x)exp(— AS") and 
G{U) = Q{AS — x) we can omit the problematic part of the AS distribution (the values 
smaller than x) from the measurement, and we get: 

(^Q,Q+iexp(-A5))g^>- 
Zq+i/Zq = , ,A5>x • (3-6) 

The price of this choice of F, G is that we do not make use of the AS" < x part of our data 
set. The value of x can be tuned to minimize the statistical error. 

Let us note that eq. (^) can be viewed as a detailed balance condition on a given U 
configuration between Q and Q + 1 sector [F and G are just the "transition probabilities"). 
This can give us a hint, that the Metropolis-step is a good a solution for F,G: F = 
min(l, exp(— AS")) and G = min(l, exp(AS')). The ratio of sectors is simply: 

y ly _ ('5o,Q+imin(l,exp(-AS')))Q 

- (<^Q,Q+imin(l,exp(A5)))g+r ^^''^ 

The inconvenient part of the distribution (AS" < 0) is cut off, however in contrast to eq. 
(|3.6| ) all configurations are used to get the expectation values. 

We have achieved our main goal: without making expensive topological sector changes 



we can obtain the ratio of sectors (see eq. (|3.5| , p^ , |3.7D ). The key point is to make 
simulations constrained to fixed topological charge, and match the results on the com- 
mon boundaries of the sectors. In the next subsection we will show in the framework of 
HMC, how to measure an expectation value, which contains a Dirac-delta on the surface. 
Our proposal is that the generation of trajectories inside a sector can be done using the 
pseudofermionic action, we do not need there the exact action. Since no sector changing is 
required, the inconvenient distribution of the pseudofermionic action jump on the bound- 
ary will not effect the measurement of the ratios of sectors. The exact action is needed 



only on the boundary: the formulas (|3.5| , p^ , |3.7|) require AS 



It is clear that using the exact fermion action when measuring the ratio of sectors 
outperforms the conventional HMC, where topological sector changing is hindered by the 
distribution of the pseudofermionic action jump. Even in that (at the moment theoretical) 
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case if we were able to use the exact fermionic action in simulations, the above presented 
method is better in determining the ratios of topological sectors. Consider a small quark 
mass simulation, where a HMC using the exact fermionic action sticks into the trivial topo- 
logical sector (now due to the fact, that nontrivial topologies are suppressed by the fermion 
determinant). If the simulation time is not long enough, then we have no information at all 
on the small (but nonvanishing) ratio Zi/Zq. However in our method this small quantity 
can be measured. The more we hit the wall from the two sides the more precisely we can 
measure Zi/Zq. The same argument applies to an R-type algorithm (where it is possible 
to use the exact jump, when crossing topological sectors [15]). 

Obviously an important issue for this new method is whether topological sectors defined 
by the overlap charge are path-connected or not. We refer to some results in the Abelian and 
in the non- Abelian gauge group case [18, 19]. If configurations with the same Q would not 
be continuously connectable in sector Q, then our assumption that we make measurements 
on the common boundary of sectors could be violated. It could happen, that the wall 
sampled from sector Q does not coincide with the wall sampled from Q -|- 1. Moreover the 
fixed sector simulations would also violate ergodicity in this case. Let us note here that the 
large autocorrelation time for the topological charge in the conventional pseudofermionic 
HMC effectively also causes the breakdown of ergodicity. In case of non-connected sectors 
one can cure these problems by releasing the system from a sector after a certain amount 
of time and closing it to another. 

Expectation value of a 5 -function 

In a HMC simulation one determines the expectation value of an operator by calcu- 
lating a sum over the measured values of the observable on N configurations, which are 
generated with the proper weights, by which they occur in the functional integral. In prac- 
tice it is not possible to measure an operator containing a Dirac-delta on the boundary 
surface on these configurations, because none of them will be exactly located on it. There- 
fore we formulate a somewhat different measurement method, and discuss its properties. 
As a result one has to make measurements at those points of the trajectories, where they 
hit the boundary. 

Let us see the details. Consider for a moment that we are able to integrate the Hamil- 
tonian equations of motion exactly. If the distribution of the gauge field, pseudofermion 
field and momenta was correct at the starting point of the trajectory, then it will be still 
correct for any of the inner points. This fact follows from energy and area conservation 
and reversibility. So we make no mistake if we put the inner points of the trajectory into 
the ensemble. We can write: 

N 1 ^ r 1 ^ 

(O) = lim — y" Oj = lim — V / 0(Ut)dt = lim — V Oj j+i (3.8) 

where t is the microcanonical time, Ut is the gauge configuration at time t, the i subscript 
at the integral means that the integration goes for the ith trajectory, and Oj^j+i is thus 
the average of the operator O along this trajectory. In the case of an observable, which 
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contains a (^-function, like those in eq. ( |3.6D or in eq. ( p.7| ) this reads 



O. 



i,i+l 



dt o{Ut) 6{\o{Ut)), 



(3.9) 



where it is indicated, that the 5-function depends on the gauge configuration only through 
the smallest absolute value eigenvalue Aq of Hw. If the variable of integration is changed 
from i to Ao, we get: 



/ 



d\[] 



dt 



dXn 



o{UtMXo). 



(3.10) 



We can go further by determining the time derivative of the smallest eigenvalue: 

IT = 9f7^ 



(3.11) 



where, again |0) is the eigenvector belonging to Aq. The trace and transpose operations 
are meant to be in color space. Recognizing, that with our previous notation 



dU 
'dt 



pU and 



(0|c/|pi|0) 



n 



yields 



dXo 
dt 



{n,p), 



dt 



dX(] 



\{n,p)\' 



N 



{O) 



lim —'S^'S^o(Ut)-r, rr 



(3.12) 



(3.13) 



tj,Xo{Uti)=0 



If we put it simple, the above formula says, that since the integration is in microcanonical 
time, the angle and velocity by which the trajectory hits the boundary has to be taken 
into account. 

Let us turn back to the case, when the integration of the equations of motion can 
be done only with finite step size integrator. The leap-frog procedure has 0(e'^) error 
per microcanonical step, which after 1/e steps makes the trajectory differ by O(e^) from 
the exact trajectory. If we can guarantee, that the modification of the trajectory at the 
boundary also violates the equations of motions only upto O(e^), then in the final results 
the errors will be proportional with 0(\/iVe^). The original reflection algorithm in [10] 
has 0(e) errors, later it was improved to O(e^) in [13]. In the Appendix we propose a 
different reflection procedure and prove that it is reversible, area-preserving and conserves 
the energy upto O(e^). 



4. Numerical results 

In the previous section we described a new method, to solve the topological sector changing 
problem of pseudofermionic HMC simulation. We describe here the details of the simula- 
tions, and finally give the topological susceptibility in physical units measured on 8^ and 
8^ X 16 lattices. 
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Simulations were done using unit length trajectories, separated by momentum and 
pseudofermion refreshments. The system was confined to a fixed topological sector in each 
run, we reflected the trajectories whenever they reached a sector boundary. The end points 
of the trajectories obviously follow the exact distribution in a given sector, usual quantities 
can be measured on them. When calculating the ratio of sectors using eq. ( |3.5| ) or eq. 
( |3.6| ) or eq. ( |3.71 ) we integrated along the trajectories, this quantity will be burdened by a 
step size error. 




Figure 2: Left panel: a typical optimization procedure of the lower limit (x) on AS in the 
formula ( |3.6| ). The statistical error of the ratio Zi/Zq shows a minimum as the function of x, which 
is considered as the optimal value. Right panel: Bare mass dependence of topological susceptibility 
using three different methods on G'' lattices. The points corresponding to the same mass were 
slightly shifted vertically for clarity. Result based on our new technique and eq. ( |3.6D is on the 
left, based on eq. ( |3.7[ ) is in the middle, the standard pseudofermionic HMC is on the right. The 
simulation parameters are from [20] 



In case of large enough statistics the value of Zq+i/Zq should be the same, indepen- 



dently which of the three formula (|3.5| , p.q , |3.7D was used to calculate it. We omit eq. 
in the following, since it is hard to give a reliable error estimate on the expectation value of 
exp(— A5), if AS can be arbitrary negative number. Eq. (|3.6D still measures exp(— AS"), 
but with a lower limit (x) on AS. Smaller limit yields a smaller and more reliable error, 
however the statistics is decreased at the same time. One can tune the value of x, so that 
the statistical error takes its minimum. A result of a typical optimum search can be seen 
on the left panel of Fig. |2|. The optimal value can be compared to the one obtained from 
eq. (|3.7|). On the right panel of Fig. || the two new topological susceptibilities and the 
one calculated by using traditional pseudofermionic HMC [20] are shown. The agreement 
is perfect. 
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m 








0.03 


3.52(13) 


0.29(11) 


2.4 


0.1 


3.17(5) 


0.53(4) 


4.3 


0.2 


2.89(2) 


0.74(6) 


5.9 


0.3 


2.88(6) 


0.99(8) 


7.9 



Table 1: Sommer-scale, pion mass and 
pion mass times box size on 8'^ x 16 lattices. 



To measure the topological susceptibility on 
8^ lattices we generated configurations with tree- 
level Symanzik improved gauge action {(3 = 4.15 
gauge coupling) and 2 step stout smeared overlap 
kernel (p = 0.15 smearing parameter, the kernel 
was the standard Wilson matrix with mo = 1.3). 
We performed runs in sectors Q = . . . 3 (based 
on the measured Z^/Z2 we can conclude, that the 
contribution oi Q > A sectors are small compared 

to statistical uncertainties). For the negatively charged sectors we used the Q — > —Q 
symmetry of the partition function. The bare masses were m = 0.03,0.1,0.2 and 0.3, 
at each mass approximately 1000 trajectories were collected. The average number of the 
topological sector boundary hittings was around 1.5 per trajectory. We calculated the ratio 
of sectors using eq. (3^) and eq. ( |3.7D . The result for the topological susceptibility can be 
seen on Fig. ^. It is nicely suppressed for the smallest mass. To convert it into physical 
units, we measured the static potential and the pion mass on 8'^ x 16 lattices (Table |l]). 
Since our statistics was quite small on these asymmetric lattices, the errors are large. Note, 
that in order to get the mass-dimension 4 topological susceptibility in physical units, one 
has to make very precise scale measurements. 




Figure 3: Topological susceptibility as the fmiction of quark mass on 8* lattices in lattice units 
(left), and in physical units as the function of pion mass (right). Scale fixing and mass measurements 
were done on 8"^ x 16. The error bars on the right plot do not contain the errors of scale fixing. The 
line is the leading order chiral behavior in the continuum. 

When interpreting the results, one should keep in mind, that the volume is small, and 
the lattice spacing is large. Note however, that smeared kernel overlap actions show nice 
scaling behavior and good locality properties already at moderate lattice spacings [21,22]. 



-10- 



5. Conclusions 

In this paper we studied the problem of topological sector changing in dynamical overlap 
simulations. The origin of the unexpectedly large autocorrelation time for the topological 
charge is connected to pseudofermions, which approximate the fermion determinant. The 
pseudofermionic action overestimates the size of the discontinuity in the fermion deter- 
minant at the topological sector boundary, so the system cannot enter easily to a new 
topological sector. This happens even if the use of the exact determinant favored a transi- 
tion. (The discontinuity of the exact determinant can be calculated in a rather inexpensive 
way.) 

We developed a new method, which circumvents the problem of topological sector 
changing. It confines the system to fixed topological sectors (by always reflecting the HMC 
trajectories from the topological sector boundaries). Thus overestimating the discontinu- 
ity of the determinant due to pseudofermions will not effect the determination of topology 
related quantities. The relative weight of two topological sector is obtained by measuring 
appropriate operators on the common boundary surface. The measurement of such oper- 
ators can be carried out by extending the usual HMC measurement method, however an 
O(e^) extrapolation in HMC step size is required. 

The new method was tested on 6'^ lattices, where previous conventional HMC results 
were available. The old and new results were consistent. We also measured the topological 
susceptibility on 8^ lattices with an improved overlap fermion and gauge action, furthermore 
simulations were done on 8^ x 16 lattices to convert the lattice results into physical units. 
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Appendix 

We present a reflection algorithm which conserves the energy upto O(e^), and which is 
reversible and area conserving. Until the boundary the trajectory is evolved by the usual 
leapfrog procedure. An elementary leapfrog step can be written in the following symbolic 
way: 



where U{e) and P(e) are the operators updating the links (u) and the momenta (p): 



where in the force term the A operator projects onto traceless, antihermitian matrices (in 
color indices). Now we split the evolution of the links into two parts: 



P(e/2)[/(e)P(e/2), 



(5.1) 



U{e) : exjp{ep)u, 



P{e) :p^p-eA u— {Sg + Spf ) 



P(e/2)C/(e/2) • C/(e/2)P(e/2). 



(5.2) 
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Consider that the boundary would be crossed during one of the evolutions of the links in 
eq. ( 15.21) . Then replace the original leapfrog with the following: 

P(e/2)[/(e/2) • U{e,)P{ec)RP{ec)Uiec) ■ U{e/2)P{e/2), 



where R is simply the reflection of the momenta (eq. (2^)). Ec is the time to reach the 
boundary surface measured from the midpoint of the leapfrog. Thus if the crossing would 
happen in the first evolution then €c < 0, if in the second, then €c > 0. The reversibility can 
be checked easily, for the energy conservation we note that both U{e)P{e) and P{e)U{e) 
conserves the energy upto 0{e^). If the evolution time does not depend on the variables, 
then both U (e) and -P(e) are area conserving. The only problem arises due to the link and 
momenta dependence of €c- We will now show that the combined effect of the €c updates 

U{ec)P{ec)RP{ec)U{ec) (5.3) 

is still area preserving. 

For simplicity we will not carry out here the computation using the SU (3) structure, 
but only for a fc-dimensional Euclidean vectorspace (the q coordinates are fc-component 
vectors). All features of the proof are also present in this simpler case. We have to 
calculate the Jacobian of the following transformation 

q'{q,p) = q + ecP + ecp' p'{q,p) = Rp + h. 

Rab = ^ab — ^HaUb is the reflection operator {a,b = 1 . . . k), where is the normal vector 
of the surface, h can be expressed as h = —CcF — ecRF, where F is the force on the wall 
(interpreted in the sector into which we reflect back the trajectory). Since h is measured 
on the wall, its q and p derivatives satisfy dh/dp = ecdh/dq (see [10]). Furthermore h is 
orthogonal to n ((n/i) = 0). For the q derivative of p' we introduce the matrix Zab^ which 
contains the derivatives of the normal vector and h: 

^P'a _ 7 

-7\ — — ^ab- 

oqb 

Using that n and h are functions of the wall coordinates, €cZ will appear in the p derivative 
of p': 

dp' 

TT-^ = Rab + ^cZab- 
OPb 

The partial derivatives of q' are: 

dqb dqb dqb 

^ = eA. + (p.+p'j|^ + e.M. 
dpb opb dpb 

With the help of the following identities (see [10]): 

dec _ ng dec _ _ rig 

dqg {np) dpa " {npY 
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the Jacobian can be written in a hypermatrix form {x o y has components XaVb)'- 



dq 



OR (np) loo/ LVi- / J 



Q is the projector to the subspace orthogonal to n (Q = 1 — n o n). J can be written as a 
product of two matrices J = J1J2, where 



J 



1 -2ecQ 
R , 



having determinant —1 (since is a reflection). 

which has an upper triangular form in the 2x2 space in a suitable basis. The determinant 
of J2 will be unity due to the orthogonality oi p + p' = p + Rp + h = —2Qp + h and n. 
Therefore det J = — 1, which shows that the transformation (^]^) preserves the area. 
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